#include <sys/types.h>
#include <sys/wait.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>

#define FIXED_POINT_SCALE_USER 1000
#define FIXED_POINT_SCALE_REAL (FIXED_POINT_SCALE_USER * FIXED_POINT_SCALE_USER)

static int fixed_point_parse(const char * input,
			     long * number)
{
	const char * parse;
	int negate;
	long value, total;
	
	parse = input;
	negate = 0;
	if (*parse == '-') {
		negate = 1;
		parse++;
	}
	
	value = strtol(parse, (char **)&parse, 10);
	if (*parse != '\0' && *parse != '.') {
		printf("%s: whole part is not a number\n", input);
		return 1;
	}
	if (*parse != '.') {
		printf("%s: no fractional component\n", input);
		return 1;
	}
	if (value < 0) {
		printf("%s: whole part is negative twice\n", input);
		return 1;
	}
	if (value > LONG_MAX / FIXED_POINT_SCALE_REAL) {
		printf("%s: whole part value is too large\n", input);
		return 1;
	}
	
	if (negate)
		value = -value;
	total = value * FIXED_POINT_SCALE_REAL;
	parse++;
	
	value = strtol(parse, (char **)&parse, 10);
	if (*parse != '\0') {
		printf("%s: fractional part is not a number\n", input);
		return 1;
	}
	if (parse < 0) {
		printf("%s: fractional part is negative\n", input);
		return 1;
	}
	if (value >= FIXED_POINT_SCALE_USER) {
		printf("%s: fractional part is too precise\n", input);
		return 1;
	}
	
	if (value != 0)
		for (; value < FIXED_POINT_SCALE_REAL / 10; value *= 10);
	if (negate)
		value = -value;
	total += value;
	
	*number = total;
	return 0;
}

static void fixed_point_print(FILE * stream,
			      long number)
{
	const char * negate;
	long whole, fract;
	
	negate = "";
	if (number < 0) {
		negate = "-";
		number = -number;
	}
	
	whole = number / FIXED_POINT_SCALE_REAL;
	fract = number % FIXED_POINT_SCALE_REAL;
	
	fprintf(stream, "%s%ld.%ld", negate, whole, fract);
}

static long fixed_point_mul(long number1,
			    long number2)
{
	return ((__int128)number1 *
		(__int128)number2) /
	       (__int128)FIXED_POINT_SCALE_REAL;
}

static long fixed_point_div(long number1,
			    long number2)
{
	return ((__int128)FIXED_POINT_SCALE_REAL *
		(__int128)number1) /
	       (__int128)number2;
}

struct coord {
	long x;
	long y;
};

static void coord_print(FILE * stream,
			const struct coord * coord)
{
	fixed_point_print(stream, coord->x);
	fprintf(stream, " ");
	fixed_point_print(stream, coord->y);
	fprintf(stream, "\n");
}

static void coord_rotate90(struct coord * coord)
{
	long xold, yold, xnew, ynew;
	
	xold = coord->x;
	yold = coord->y;
	
	xnew = -yold;
	ynew = xold;
	
	coord->x = xnew;
	coord->y = ynew;
}

static void coord_rotate180(struct coord * coord)
{
	long xtemp, ytemp;
	
	xtemp = coord->x;
	ytemp = coord->y;
	
	xtemp = -xtemp;
	ytemp = -ytemp;
	
	coord->x = xtemp;
	coord->y = ytemp;
}

static int coord_nonzero(const struct coord * coord)
{
	return (coord->x | coord->y) != 0;
}

static void coord_sub(const struct coord * coord1,
		      const struct coord * coord2,
		      struct coord * result)
{
	result->x = coord1->x - coord2->x;
	result->y = coord1->y - coord2->y;
}

static void coord_mid(const struct coord * coord1,
		      const struct coord * coord2,
		      struct coord * result)
{
	result->x = (coord1->x + coord2->x) >> 1;
	result->y = (coord1->y + coord2->y) >> 1;
}

static int coord_inline(const struct coord * coord1,
			const struct coord * coord2)
{
	return fixed_point_mul(coord1->x, coord2->y) ==
	       fixed_point_mul(coord1->y, coord2->x);
}

struct point {
	struct coord point;
	int valid;
};

struct line {
	struct coord point;
	struct coord slope;
};

struct rectangle {
	long xmax;
	long xmin;
	long ymax;
	long ymin;
};

static void line_rotate90(struct line * line)
{
	coord_rotate90(&line->slope);
}

static void line_rotate180(struct line * line)
{
	coord_rotate180(&line->slope);
}

static void line_through(struct line * line,
			 const struct coord * coord1,
			 const struct coord * coord2)
{
	coord_mid(coord1, coord2, &line->point);
	coord_sub(coord1, coord2, &line->slope);
}

static int line_contains(const struct line * line,
			 const struct coord * coord)
{
	struct coord delta;
	
	coord_sub(&line->point, coord, &delta);
	
	return coord_inline(&line->slope, &delta);
}

static int line_equal(const struct line * line1,
		      const struct line * line2)
{
	struct coord delta;
	
	coord_sub(&line1->point, &line2->point, &delta);
	if (coord_nonzero(&delta)) {
		if (!coord_inline(&delta, &line1->slope))
			return 0;
		if (!coord_inline(&delta, &line2->slope))
			return 0;
	} else {
		if (!coord_inline(&line1->slope, &line2->slope))
			return 0;
	}
	
	return 1;
}

static void line_escapes(const struct line * line,
			 const struct rectangle * bounds,
			 struct coord * coord)
{
	long xbound, ybound, xticks, yticks;
	
	xbound = 0;
	ybound = 0;
	xticks = LONG_MAX;
	yticks = LONG_MAX;
	
	if (line->slope.x != 0) {
		if (line->slope.x >= 0)
			xbound = bounds->xmax;
		else
			xbound = bounds->xmin;
		xticks = fixed_point_div(xbound - line->point.x, line->slope.x);
	}
	if (line->slope.y != 0) {
		if (line->slope.y >= 0)
			ybound = bounds->ymax;
		else
			ybound = bounds->ymin;
		yticks = fixed_point_div(ybound - line->point.y, line->slope.y);
	}
	
	if (xticks < yticks)
		ybound = fixed_point_mul(xticks, line->slope.y) + line->point.y;
	if (yticks < xticks)
		xbound = fixed_point_mul(yticks, line->slope.x) + line->point.x;
	
	coord->x = xbound;
	coord->y = ybound;
}

struct array_of_points {
	struct point * array;
	unsigned int count;
};

struct array_of_lines {
	struct line * array;
	unsigned int count;
};

static int initialize_data(struct array_of_points * points,
			   struct array_of_lines * lines,
			   int argc, char ** argv)
{
	unsigned int count, i;
	
	if ((argc & 1) == 0 ||
	    argc <= 1) {
		printf("%s <<point x> <point y>> ...\n", argv[0]);
		return 1;
	}
	
	argc--;
	argv++;
	
	count = argc >> 1;
	points->array = malloc(count * sizeof(*points->array));
	for (i = 0; i < count; i++) {
		struct coord * coord;
		
		coord = &points->array[i].point;
		if (fixed_point_parse(*argv, &coord->x))
			return 1;
		argv++;
		if (fixed_point_parse(*argv, &coord->y))
			return 1;
		argv++;
	}
	points->count = count;
	
	count *= count - 1;
	lines->array = malloc(count * sizeof(*lines->array));
	lines->count = 0;
	
	return 0;
}

static void append_candidate(struct array_of_lines * lines,
			     const struct line * candidate)
{
	unsigned int i;
	
	for (i = 0; i < lines->count; i++)
		if (line_equal(lines->array + i, candidate))
			return;
	
	lines->array[i] = *candidate;
	lines->count++;
}

static void find_symmetry(struct array_of_points * points,
			  struct array_of_lines * lines)
{
	unsigned int i, j, k, last;
	
	for (i = 0; i < points->count - 1; i++) {
		for (j = i + 1; j < points->count; j++) {
			const struct coord * coord1;
			const struct coord * coord2;
			struct line line;
			
			coord1 = &points->array[i].point;
			coord2 = &points->array[j].point;
			
			line_through(&line, coord1, coord2);
			append_candidate(lines, &line);
			line_rotate90(&line);
			append_candidate(lines, &line);
		}
	}
	
	last = 0;
	
	for (k = 0; k < lines->count; k++) {
		const struct line * line;
		int valid;
		
		line = lines->array + k;
		for (i = 0; i < points->count; i++)
			points->array[i].valid = 0;
		valid = 1;
		
		for (i = 0; i < points->count; i++) {
			const struct coord * coord1;
			int buddy;
			
			if (points->array[i].valid)
				continue;
			coord1 = &points->array[i].point;
			if (line_contains(line, coord1))
				continue;
			buddy = 0;
			
			for (j = i + 1; j < points->count; j++) {
				const struct coord * coord2;
				struct line cut;
				
				if (points->array[j].valid)
					continue;
				coord2 = &points->array[j].point;
				line_through(&cut, coord1, coord2);
				line_rotate90(&cut);
				
				if (line_equal(line, &cut)) {
					points->array[j].valid = 1;
					buddy = 1;
					break;
				}
			}
			
			if (!buddy) {
				valid = 0;
				break;
			}
		}
		
		if (valid) {
			if (last != k)
				lines->array[last] = *line;
			last++;
		}
	}
	
	lines->count = last;
}

static char filename_ldata[] = "ldata-XXXXXX";
static char filename_pdata[] = "pdata-XXXXXX";
static char filename_gpcmd[] = "gpcmd-XXXXXX";

static FILE * create_tempfile(char * template)
{
	int fd;
	FILE * fp;
	
	fd = mkstemp(template);
	fp = fdopen(fd, "w");
	
	return fp;
}

static void output_results(const struct array_of_points * points,
			   const struct array_of_lines * lines)
{
	FILE * ldata, * pdata, * gpcmd;
	struct rectangle bounds;
	unsigned int i;
	long xrange, yrange, trange;
	pid_t pid;
	
	ldata = create_tempfile(filename_ldata);
	pdata = create_tempfile(filename_pdata);
	gpcmd = create_tempfile(filename_gpcmd);
	
	bounds.xmax = LONG_MIN;
	bounds.xmin = LONG_MAX;
	bounds.ymax = LONG_MIN;
	bounds.ymin = LONG_MAX;
	for (i = 0; i < points->count; i++) {
		const struct coord * coord;
		
		coord = &points->array[i].point;
		if (bounds.xmax < coord->x)
			bounds.xmax = coord->x;
		if (bounds.xmin > coord->x)
			bounds.xmin = coord->x;
		if (bounds.ymax < coord->y)
			bounds.ymax = coord->y;
		if (bounds.ymin > coord->y)
			bounds.ymin = coord->y;
	}
	
	xrange = bounds.xmax - bounds.xmin;
	yrange = bounds.ymax - bounds.ymin;
	trange = 0;
	if (trange < xrange)
		trange = xrange;
	if (trange < yrange)
		trange = yrange;
	trange <<= 1;
	xrange = trange - xrange;
	yrange = trange - yrange;
	xrange >>= 1;
	yrange >>= 1;
	bounds.xmax += xrange;
	bounds.xmin -= xrange;
	bounds.ymax += yrange;
	bounds.ymin -= yrange;
	
	for (i = 0; i < lines->count; i++) {
		struct line line;
		struct coord coord;
		
		if (i != 0)
			fprintf(ldata, "\n");
		
		line = lines->array[i];
		line_escapes(&line, &bounds, &coord);
		coord_print(ldata, &coord);
		line_rotate180(&line);
		line_escapes(&line, &bounds, &coord);
		coord_print(ldata, &coord);
	}
	
	for (i = 0; i < points->count; i++) {
		const struct coord * coord;
		
		if (i != 0)
			fprintf(pdata, "\n");
		
		coord = &points->array[i].point;
		coord_print(pdata, coord);
	}
	
	fprintf(gpcmd, "set terminal wxt size 800,800\n");
	fprintf(gpcmd, "set xrange [");
	fixed_point_print(gpcmd, bounds.xmin);
	fprintf(gpcmd, ":");
	fixed_point_print(gpcmd, bounds.xmax);
	fprintf(gpcmd, "]\n");
	fprintf(gpcmd, "set yrange [");
	fixed_point_print(gpcmd, bounds.ymin);
	fprintf(gpcmd, ":");
	fixed_point_print(gpcmd, bounds.ymax);
	fprintf(gpcmd, "]\n");
	fprintf(gpcmd, "unset key\n");
	fprintf(gpcmd, "plot '%s' with lines, '%s' with points\n",
		filename_ldata, filename_pdata);
	
	fclose(ldata);
	fclose(pdata);
	fclose(gpcmd);
	
	pid = fork();
	if (pid == 0) {
		execlp("gnuplot",
		       "gnuplot",
		       "-p",
		       filename_gpcmd,
		       (char *)NULL);
		printf("gnuplot: failed to execute\n");
		exit(-1);
	} else {
		waitpid(pid, NULL, 0);
	}
	
	unlink(filename_ldata);
	unlink(filename_pdata);
	unlink(filename_gpcmd);
}

int main(int argc, char ** argv)
{
	struct array_of_points points;
	struct array_of_lines lines;
	
	if (initialize_data(&points, &lines, argc, argv))
		return 0;
	find_symmetry(&points, &lines);
	output_results(&points, &lines);
	
	return 0;
}
